using System;
using System.Diagnostics;
using System.Globalization;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;

namespace Atomic.Structures
{
	[Serializable]
	[DebuggerDisplay("{DebuggerDisplay}")]
	public class Basis
	{
		[NonSerialized] int hashCode;

		static Basis()
		{
			// Special constructor for the standard basis.
			StandardBasis = new Basis();
		}

		/// <summary>
		/// Creates a new basis with basis vectors specified relative to the standard basis (i.e. physical coordinates).
		/// </summary>
		public Basis(SpaceVector a, SpaceVector b, SpaceVector c)
			: this(new Position(a), new Position(b), new Position(c))
		{
		}

		/// <summary>
		/// Creates a new basis with basis vectors specified relative to a given basis.
		/// </summary>
		public Basis(Basis parentBasis, SpaceVector a, SpaceVector b, SpaceVector c)
			: this(new Position(parentBasis, a), new Position(parentBasis, b), new Position(parentBasis, c))
		{
		}

		public Basis(Position a, Position b, Position c)
		{
			if (a == null || b == null || c == null)
			{
				throw new ArgumentNullException();
			}

			Basis basis = a.Basis;
			if (b.Basis != basis || c.Basis != basis)
			{
				// No common basis. Transform to the standard basis.
				basis = Basis.StandardBasis;
				a = new Position(a.PhysicalCoordinates);
				b = new Position(b.PhysicalCoordinates);
				c = new Position(c.PhysicalCoordinates);
			}

			ParentBasis = basis;
			A = a;
			B = b;
			C = c;

			SpaceMatrix map = new SpaceMatrix(a.Coordinates, b.Coordinates, c.Coordinates);
			SpaceMatrix inverseMap;
			if (!SpaceMatrix.TryInverse(map, out inverseMap))
			{
				throw new ArgumentException("Vectors aren't linearly independent.");
			}
			Map = map;
			InverseMap = inverseMap;

			if (basis == Basis.StandardBasis)
			{
				// Don't perform arithmetics to avoid unnecessary round-off errors.
				PhysicalMap = map;
				InversePhysicalMap = inverseMap;
			}
			else
			{
				PhysicalMap = basis.PhysicalMap * map;
				InversePhysicalMap = inverseMap * basis.InversePhysicalMap;
			}
		}

		/// <summary>
		/// Creates a new basis with lengths of basis vectors and Euler angles given. Angles are in degrees following the crystallographic convention.
		/// </summary>
		public Basis(double a, double b, double c, double alpha, double beta, double gamma)
			: this(new SpaceVector(a, 0.0, 0.0),
				new SpaceVector(b * Math.Cos(gamma * Math.PI / 180), b * Math.Sin(gamma * Math.PI / 180), 0.0),
				new SpaceVector(c * Math.Cos(beta * Math.PI / 180), 0.0, c))
		{
			// http://en.wikipedia.org/wiki/Basis_%28linear_algebra%29#Expression_of_a_basis

			// self._vecc = [0.,0.,0.]
			// self._vecc[0] = self._c * cos( betar )
			// cy = ( cos(alphar) - cos(gammar) * cos(betar) ) / sin(gammar)
			// self._vecc[1] = self._c * cy
			// cz = sqrt( ( sin( betar ) )**2 - cy**2 )
			// self._vecc[2] = self._c * cz

			if (a <= 0.0 || b <= 0.0 || c <= 0.0)
			{
				throw new ArgumentException();
			}

			if ((alpha != 90.0 || beta != 90.0 || gamma != 90.0) && (alpha != 90.0 || beta != 90.0 || gamma != 120.0 || a != b))
			{
				throw new NotImplementedException();
			}

			//if (alpha < 0.0 || alpha > 180.0 || ... throw new ArgumentException()
		}

		private Basis()
		{
			// Special constructor for the standard basis.
			Map = SpaceMatrix.Identity;
			InverseMap = SpaceMatrix.Identity;
			PhysicalMap = SpaceMatrix.Identity;
			InversePhysicalMap = SpaceMatrix.Identity;
			ParentBasis = this;
			A = new Position(this, 1.0, 0.0, 0.0);
			B = new Position(this, 0.0, 1.0, 0.0);
			C = new Position(this, 0.0, 0.0, 1.0);
		}

		public static Basis CommonBasis(Basis basis1, Basis basis2, Basis basis3)
		{
			if (basis1 == null || basis2 == null || basis3 == null)
			{
				throw new ArgumentNullException();
			}

			Basis basis = basis1;
			if (basis2 != basis || basis3 != basis)
			{
				// No common basis.
				basis = Basis.StandardBasis;
			}

			return basis;
		}

		/// <summary>
		/// Writes the position in coordinates of this basis.
		/// </summary>
		public Position Decompose(Position position)
		{
			if (position.Basis == this)
			{
				// Already using the right basis.
				return position;
			}

			if (this == Basis.StandardBasis)
			{
				// Don't perform arithmetics to avoid unnecessary round-off errors.
				return new Position(position.PhysicalCoordinates);
			}
			else
			{
				return new Position(this, InversePhysicalMap * position.PhysicalCoordinates);
			}
		}

		/// <summary>
		/// Writes the position in coordinates of this basis. Rounds off integer values.
		/// </summary>
		public Position DecomposeRounded(Position position)
		{
			double tol = 1.0e-10;
			Position p = Decompose(position);

			double x = p.Coordinates.X;
			if (Math.Abs(x - Math.Round(x)) < tol)
			{
				x = (int)Math.Round(x);
			}

			double y = p.Coordinates.Y;
			if (Math.Abs(y - Math.Round(y)) < tol)
			{
				y = (int)Math.Round(y);
			}

			double z = p.Coordinates.Z;
			if (Math.Abs(z - Math.Round(z)) < tol)
			{
				z = (int)Math.Round(z);
			}

			return new Position(this, x, y, z);
		}

		public Basis Scale(double scaling)
		{
			if (scaling <= 0.0)
			{
				throw new ArgumentException();
			}

			return new Basis(A * scaling, B * scaling, C * scaling);
		}

		public Basis Rotate(SpaceMatrix rotation)
		{
			if (Math.Abs(SpaceMatrix.Determinant(rotation) - 1.0) > 1.0e-6)
			{
				throw new ArithmeticException("Invalid rotation matrix.");
			}

			return new Basis(rotation * A.PhysicalCoordinates, rotation * B.PhysicalCoordinates, rotation * C.PhysicalCoordinates);
		}

		public virtual bool Equals(Basis other)
		{
			if (object.ReferenceEquals(other, null))
			{
				return false;
			}

			//return EqualityHelper.Equals(this, other, new object[] { ParentBasis, A, B, C }, new object[] { other.ParentBasis, other.A, other.B, other.C }, hashCode, other.hashCode);
			return EqualityHelper.Equals(this, other, new object[] { PhysicalMap }, new object[] { other.PhysicalMap }, hashCode, other.hashCode);
		}

		public override bool Equals(object other)
		{
			return Equals(other as Basis);
		}

		public override int GetHashCode()
		{
			if (hashCode == 0)
			{
				//hashCode = EqualityHelper.HashCode(new object[] { ParentBasis, A, B, C });
				hashCode = EqualityHelper.HashCode(new object[] { PhysicalMap });
			}

			return hashCode;
		}
		
		public static bool operator ==(Basis basis1, Basis basis2)
		{
			if (object.ReferenceEquals(basis1, null))
			{
				if (object.ReferenceEquals(basis2, null))
				{
					return true;
				}

				return false;
			}

			return basis1.Equals(basis2);
		}

		public static bool operator !=(Basis basis1, Basis basis2)
		{
			return !(basis1 == basis2);
		}

		/// <summary>
		/// The basis vectors themselves could be specified relative to a common parent basis as given here.
		/// If no parent basis is specified, this is just the standard basis.
		/// </summary>
		public Basis ParentBasis
		{
			get;
			private set;
		}

		/// <summary>
		/// First basis vector.
		/// </summary>
		public Position A
		{
			get;
			private set;
		}
		
		/// <summary>
		/// Second basis vector.
		/// </summary>
		public Position B
		{
			get;
			private set;
		}

		/// <summary>
		/// Third basis vector.
		/// </summary>
		public Position C
		{
			get;
			private set;
		}

		/// <summary>
		/// Matrix defining the mapping of a vector relative to this basis to the parent basis.
		/// </summary>
		public SpaceMatrix Map
		{
			get;
			private set;
		}

		public SpaceMatrix InverseMap
		{
			get;
			private set;
		}

		/// <summary>
		/// Matrix defining the mapping of a vector relative to this basis to the standard basis (i.e. physical coordinates).
		/// </summary>
		public SpaceMatrix PhysicalMap
		{
			get;
			private set;
		}

		public SpaceMatrix InversePhysicalMap
		{
			get;
			private set;
		}

		public static Basis StandardBasis
		{
			get;
			private set;
		}

		private string DebuggerDisplay
		{
			get
			{
				return string.Format(CultureInfo.InvariantCulture, "A = ({0:f4}, {1:f4}, {2:f4}), B = ({3:f4}, {4:f4}, {5:f4}), C = ({6:f4}, {7:f4}, {8:f4})",
					A.PhysicalCoordinates.X, A.PhysicalCoordinates.Y, A.PhysicalCoordinates.Z, 
					B.PhysicalCoordinates.X, B.PhysicalCoordinates.Y, B.PhysicalCoordinates.Z, 
					C.PhysicalCoordinates.X, C.PhysicalCoordinates.Y, C.PhysicalCoordinates.Z);
			}
		}
	}
}
